function [stored_hist stored_diams ThresholdRange diams_range,area_range] = ThresholdingProcess(depth, NoDilate, ThresholdRange,diams_range,area_range)
%% ThresholdingProcess perform segmentation and glomeruli counting at various threshold levels.
%
%  Arguments:
%   depth    3D image (single or double)
%   NoDilate (optional) binary 0=no dilation or 1=run dilation
%   ThresholdRange  (optional) array of threshold values eg. 0.11:0.01:0.5
%   diams_range    (optional) 2D array of diamter histograms at each
%  threshold
%   area_range    (optional) 2D array of 'areas' or size of labels
%
%  Usage: 
%   >> load Cache/1008.2.18.3.1.1.1/depth_mask_01.mat
%   >> [stored_hist stored_diams] = ThresholdingProcess(depth, 1);
%   >> PlotXGlomOutputDist(stored_hist, stored_diams,'1008.2.18.3.1.1.1',0.0508);
%
%
%   - Michael Eager,  (michael.eager@monash.edu)
%   - Monash Biomedical Imaging,  Monash University, Australia

%     Copyright © 2012-2013 Michael Eager <michael.eager@monash.edu> 
%
%     This file is part of Xglom.
%
%     This is free software: you can redistribute it and/or modify
%     it under the terms of the GNU General Public License as published by
%     the Free Software Foundation, either version 3 of the License, or
%     (at your option) any later version.
%
%     This is distributed in the hope that it will be useful,
%     but WITHOUT ANY WARRANTY; without even the implied warranty of
%     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
%     GNU General Public License for more details.
%
%     You should have received a copy of the GNU General Public License
%     along with this program.  If not, see <http://www.gnu.org/licenses/>.

Xglompath = which('Xglom');
[Xglom_path,Xglom_filename,MRtools_ext] = fileparts(Xglompath);
addpath(fullfile(Xglom_path,'plugins','progressbar'));
tstart=tic;

% Set up progress bar with abort button
gui_active(1);      % will add an abort button
hDummy=figure;
hProgress=figure;
hProgress_0 = hProgress;
hProgress  = progressbar([],0,' Performing Threshold Statistics (50 steps)' );

close(hProgress_0);
close(hDummy)
telapsed = toc(tstart);

if nargin < 2
NoDilate= 0;
end
if nargin < 3
   ThresholdRange= 0.11:0.01:0.5;
end
if nargin < 4
diams_range=1:0.1:15;
end
if nargin < 5
    area_max = 1000; %max(areas);
    area_range=(1:area_max);
end

stored_hist=zeros(length(ThresholdRange),length(area_range));
stored_diams = zeros(length(ThresholdRange),length(diams_range));
counter=1;
incr_progress = double(1)/double(length(ThresholdRange)); 

for counter=1:length(ThresholdRange) % Loop over different Thresholds
    %% Thresholding
    Holes=depth>ThresholdRange(counter);
    
    %% Discard single voxels
    if NoDilate
     ConArray=conndef(3,'minimal');
     ConArray(2,2,2)=0;
     Holes=imdilate(Holes,ConArray).*Holes;
     Holes=single(Holes);
    end
    
    %% Statistics, Label Glomerular Regions
    CC = bwconncomp(Holes,6);

    
    %% calculate area histogram
    areas = cellfun(@numel,CC.PixelIdxList);
    [area_hist, area_axis]=hist(areas,area_range);
   
    stored_hist(counter,:) = area_hist;
    disp([' Total number of detected objects in area range: ', num2str(sum(area_hist(1:end-1)))])
    
    % Diameter histogram from area data assuming spherical glomeruli
    
    [diams ~]=hist((6/pi*areas(areas>1)).^(1/3),diams_range);
    diams(end)=0;
    stored_diams(counter,:) = diams;
    
    telapsed = toc(tstart);

    if ~hProgress
        return
    end
    if counter == length(ThresholdRange)
        
        hProgress = progressbar( hProgress, incr_progress, 'Completed Xglom thresholding process.');
    else
        
        hProgress = progressbar( hProgress,incr_progress , ['Threshold process (step ' num2str(counter) ' of ' num2str(length(ThresholdRange)) ')']);
    end
end
telapsed = toc(tstart);
disp(['Processing time (sec): ' num2str(telapsed)])

hProgress = progressbar( hProgress, -1 , 'Image Processing Complete. Processing output plots.');

